///////////////////////////////////////////////////////////////////////////////
// Copyright (C) Nial Peters 2012
//
// This file is part of AvoScan.
//
// AvoScan is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// AvoScan is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with AvoScan.  If not, see <http://www.gnu.org/licenses/>.
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


#include"Python.h"
#include "structmember.h"
#include<math.h>
#include"autodealloc_numpy.h"

//macro for converting angles in degrees to radians
#define RADIANS(x) ((x) * 0.017453292519943295)

static char mod_doc[] = "The cPlume module provides a few functions used by ica_inspect \
that were just too slow when implemented in pure Python. You probably don't want to use any \
of them yourself.";



static char plumeDepth_doc[] = "By assuming a cylindrical plume, this function \
calculates the the depth of plume being looked through for a given set of scan \
parameters. get_path_length(r,theta0,h_dist,v_angle,theta) where r is the radius \
of the plume (of the cylinder), theta0 is the scan angle of the centre of the plume, \
h_dist is the horizontal distance (in metres) of the scanner from the centre of the \
plume, v_angle is the angle of the scan plane from horizontal and theta is a 1D numpy \
array of scan angles that you want the scan path calculated for. Scan angles which do \
not intersect the plume or are tangent to the plume will return 0 as path length.";



static PyObject * cPlume_getPathLength(PyObject *self, PyObject *args){

	double r; //plume radius (radius of cylinder)
	double theta0; //angle of the centre of the plume
	double h_dist; //horizontal distance from the plume
	double v_angle; //angle of scanner
	PyObject *theta_arr; //numpy array of theta values

	double cx1,cx2,cy1,cy2;
    double c,d;
    double A,B,C, discriminant, sqrt_discr;
    double theta;
	int i=0;//counter

	//parse the arguments passed to the function by Python - no increase in
	//ref count
	if(!PyArg_ParseTuple(args, "ddddO", &r, &theta0, &h_dist, &v_angle, &theta_arr)){
		PyErr_SetString(PyExc_ValueError,"Invalid parameters");
		return NULL;
	}

	//check that we have been passed a 1D numpy array for theta values
	if (!PyArray_Check(theta_arr)){
		PyErr_SetString(PyExc_ValueError,"theta must be a numpy array");
		return NULL;
	}
	if (PyArray_NDIM(theta_arr) != 1){
		PyErr_SetString(PyExc_ValueError,"theta must only have one dimension");
		return NULL;
	}

	//get the size of the theta array
	int num_entries =(int)PyArray_DIM(theta_arr, 0);

	//ellipse: x^2/a^2 + y^2/b^2 = 1
    //ellipse params are constrained by the viewing angle
    double cos_v_angle = cos(RADIANS(v_angle));
	double a = r;
    double b = r / cos_v_angle;

    //transform h_dist into the plane of the ellipse
    h_dist = h_dist / cos_v_angle;

    //allocate memory for results array
    double *result = malloc(num_entries*sizeof(double));
    if (result == NULL){
    	PyErr_SetString(PyExc_ValueError,"Could not alloc results array");
    	return NULL;
    }

    //do the calculation!
    for(i; i<num_entries; i++){

    	theta = *((double*)PyArray_GETPTR1(theta_arr,i));

    	//use the equation of a straight line x = cy + d we do it
		//this way around because it is possible that theta==theta0, which
		//leads to x1==x2 and an infinite gradient. However, h_dist is never
		//0 and so dx/dy will never be infinite.
		c = tan(RADIANS(theta-theta0));
		d = c * h_dist;

		//substitute the expression for x from the line into the equation of the ellipse
		//and solve using the quadratic equation.
		//quadratic equation -> x = (-B +- sqrt(B^2 - 4AC))/(2A)
		A = ((a * a) / (b * b)) + (c * c);
		B = 2.0 * c * d;
		C = (d * d) - (a * a);

		//mask out scan angles that did not intersect the plume (they have imaginary roots)
		discriminant = (B * B) - (4 * A * C);
		if (discriminant <= 0.0){
			result[i] = 0;
			continue;
		}

		//solve for the intersection points of the lines with the ellipse
		sqrt_discr = sqrt(discriminant);
		cy1 = (-B + sqrt_discr) / (2.0 * A);
		cy2 = (-B - sqrt_discr) / (2.0 * A);

		//put the y values back into the equation of a straight line and solve for x
		cx1 = c * cy1 + d;
		cx2 = c * cy2 + d;

		result[i] = sqrt(((cx2 - cx1) *(cx2 - cx1))  + ((cy2 - cy1) *(cy2 - cy1)));
    }

    //create numpy array as return value
    npy_intp dims[1];
    dims[0] = num_entries;
    return createAutoDeallocNumpyArray (1, dims, NPY_FLOAT64, result);

}



//set up the functions to be visible in Python
static PyMethodDef cPlume_methods[] = {
    {"get_path_length", cPlume_getPathLength, METH_VARARGS, plumeDepth_doc},
    {NULL, NULL}
};



//set up module to be importable in Python
PyMODINIT_FUNC initcPlume(void){
    import_array();
    import_autodealloc_numpy();
    Py_InitModule3("cPlume", cPlume_methods, mod_doc);
}
